%%牛顿插值函数

%xi、yi为已知的数据
%x为要预测的数据点
function p = Newton_fun(x,xi,yi)
n=length(xi);  %获取输入数据点的个数
f=zeros(n,n);  %
q=zeros(1,n);

% 将均差表的第一列赋值yi
for k=1:n      
    f(k,1)=yi(k);
end

% 求均差表
for i=2:n  %i表示均差结束      
    for k=i:n
        %计算各阶的均差值
        f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));  
    end
end

disp('均差表如下：');
for i=1:n       
    for k=1:n
        %以保留小数点后4位有效数字的方式打印均差表
        fprintf("%2.4f  ", f(i,k)) ;
    end
    fprintf('\n');
end

%%求插值多项式各项的系数
for i = 1:1:n
    q(i)=f(i,i); %q表示牛顿插值法的各项系数
end
disp('牛顿插值多项式各项的系数为：\n');
for i = 1:1:n
    fprintf('q[%d] = %.6f  ', i, q(i));
end

nn = length(x); %预测数据的个数
p = zeros(1,nn); %p表示预测值

%%牛顿插值法预测的数据值计算
for i = 1:1:nn
    for j = 1:1:n
        factorial = 1;
        for k = 1:1:j - 1
            factorial = factorial*(x(i)-xi(k));
        end
        factorial = q(j)*factorial;
        p(i) = p(i) + factorial;
    end
end

end
